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In the present study we generaUze the self-consistent Hartree-Fock-Bogohubov (HFB) theory 
formulated in the coordinate space to the case which incorporates an arbitrary mixing between 
protons and neutrons in the particle-hole (p-h) and particle-particle (p-p or pairing) channels. We 
define the HFB density matrices, discuss their spin-isospin structure, and construct the most general 
energy density functional that is quadratic in local densities. The consequences of the local gauge 
invariance are discussed and the particular case of the Skyrme energy density functional is studied. 
By varying the total energy with respect to the density matrices the self-consistent one-body HFB 
Hamiltonian is obtained and the structure of the resulting mean fields is shown. The consequences 
of the time-reversal symmetry, charge invariance, and proton-neutron symmetry are summarized. 
The complete list of expressions required to calculate total energy is presented. 
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I. INTRODUCTION 

One of the main goals of nuclear theory is to build 
the unified microscopic framework for heavy nuclei in 
which the bulk nuclear properties, nuclear excitations, 
and nuclear reactions can be described on the same foot- 
ing. Microscopic theory also provides the solid founda- 
tion for phenomenological models and coupling schemes 
which have been applied so successfully to explain specific 
nuclear properties. Exotic short-lived nuclei are very im- 
portant in this quest. The abnormal ncutron-to- proton 
ratios of these nuclei isolate and amplify important fea- 
tures, which are not clearly visible in stable systems. 

For medium-mass and heavy nuclei, a critical challenge 
is the quest for the universal energy density functional, 
which will be able to describe properties of finite nuclei 
as well as extended asymmetric nucleonic matter (e.g., as 
found in neutron stars). Self-consistent methods based 
on the density functional theory have already achieved a 
level of sophistication and precision which allows analyses 
of experimental data for a wide range of properties and 
for arbitrarily heavy nuclei. For instance, self-consistent 
HE and HFB models are now able to reproduce measured 
nuclear binding energies with an impressive rms error of 
^700 keV [1, 2]. However, much work remains to be 
done. Developing a universal nuclear density functional 
will require a better understanding of the density depen- 
dence, isospin effects, pairing, as well as an improved 
treatment of many-body correlations. All those aspects 
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are essential for the structure of proton-rich nuclei with 
NKiZ, which are expected to exhibit proton-neutron (pn) 
pairing [3]; it is precisely in those nuclei that the state- 
of-the-art microscopic mass formula needs to be supple- 
mented by a phenomenological Wigner term [1, 2]. 

In spite of an impressive experimental progress in the 
heavy Nk,Z region, it is still unclear (i) what the spe- 
cific fingerprints of the pn pairing are and (ii) what is 
the interplay between the like-particle and pn (r=0,l) 
p-h and p-p channels. Before attempting to answer these 
questions, established theoretical models of nuclear pair- 
ing need to be generalized to properly account for pn 
correlations. The present work is a step in this direction. 
We propose the general HFB formalism which fully in- 
corporates the pn mixing on the mean-field level. The 
resulting density matrices have a very rich spin-isospin 
structure, which, in the presence of static pn pairing, 
can produce novel mean- fields and deformations. 

The paper is organized as follows. Section II contains 
a brief review of the pn pairing. Section III discusses 
the density matrices (scalar, vector, and tensor), both 
in the p-h and p-p channel. The discussion is based on 
the coordinate-space HFB formalism [4-6] , which was in- 
troduced earlier to describe pairing correlations between 
like nucleons. This method is the tool of choice when 
dealing with weakly-bound heavy nuclei [7]. The energy 
functional is constructed in Sec. IV, the associated mean 
fields are derived in Sec. V, and Sec. VI deals with the 
resulting coordinate-space HFB equations. In the discus- 
sion of pn pairing, the notion of self-consistent symme- 
tries, especially those associated with charge invariance 
and time reversal, is crucial, and we devote Sec. VII to 
this topic. Finally, conclusions are contained in Sec. VIII. 
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II. PROTON-NEUTRON PAIRING, A CONCISE 
OVERVIEW 

A unique aspect of proton-rich nuclei with NwZ is that 
neutrons and protons occupy the same single-particle or- 
bitals. Consequently, due to the large spatial overlaps be- 
tween neutron and proton single-particle wave functions, 
pn pairing is expected to be present in those systems. 

So far, the strongest evidence for enhanced pn corre- 
lations around the N=Z line comes from the measured 
binding energies [8 16] and the isospin structure of the 
low- lying states in odd-odd nuclei [16-28] . The pn corre- 
lations are also expected to play some role in single-beta 
decay [29 31], double-beta decay [32-38], transfer reac- 
tions [39 45] (sec, however, Ref. [16]), structure of low- 
lying collective states [46] , alpha decay and alpha correla- 
tions [45, 47-52], structure of high spin states [14, 20, 53- 
79], and in properties of low-density nuclear matter [80- 
90]. 

Actually, the pn pairing is not "the new kid on the 

block" but it has a long history and is ultimately con- 
nected to the charge invariancc of the strong Hamilto- 
nian. (For reference, in 1932 Heisenberg introduced iso- 
topic spin [91] and in 1936 Wigner introduced the nuclear 
SU(4) supermultiplets [92].) An important step was the 
adaptation of Racah's concept of seniority by Racah and 
Talmi [93], and Flowers [94] in 1952. In the indepen- 
dent quasiparticle (BCS) picture [95], pairing condensate 
appears as a result of an attractive interaction between 
quasiparticles near the Fermi surface. The term "nuclear 
superconductivity" was first used by Pines at the 1957 
Rehovot Conference to point out that the new BCS the- 
ory might also apply to nuclei [96]. This was formally 
accomplished in the late fifties [97, 98] and shortly af- 
terwards the importance of pn pairing was emphasized 
[47, 99, 100] and a number of theoretical papers deal- 
ing with the generalization of the BCS theory to the pn 
pairing case appeared [101-103]. 

Independently, group-theoretical methods based on the 
quasi-spin formalism were developed. Many insights were 
gained by simple solvable models employing symmetry- 
dictated interactions [104-111]. Two families of mod- 
els were used, one based on the j-j coupling with the 
symmetry SO (5) (appropriate for the T=l pairing) and 
the other based on the L-S coupling with the symmetry 
S0(8) (appropriate for the T=0 and T=l pairing). These 
models have been consecutively developed and applied to 
various physically interesting cases [36, 42, 44, 45, 112- 
114]. Among many techniques used to solve the problem 
of pn pairing with schematic interactions, worth mention- 
ing are the exact methods [43, 115, 116] used to describe 
isovector states of a charge-independent pairing Hamil- 
tonian. 

Properties of pn pairing (at low and high spins, in cold 
and hot nuclei) have been studied within the large-scale 
shell model (diagonalization shell-model, variational shell 
model, and Monte Carlo shell model) [20, 21, 64, 69, 70, 
90, 117-123]. It was concluded that the isovector pair- 



ing in the dominating J=0 channel mainly acts between 
time-reversed states within the same shell. On the other 
hand, isoscalar pairing can also involve coupling (mainly 
J=l) between spin-orbit partners. Consequently, spin- 
orbit splitting plays a crucial role in understanding the 
T=0 pairing [20]. 

It is to be noted that it is by no means obvious how 
to extract "pairing correlations" from the realistic shell- 
model calculations. The "pairing Hamiltonian" is an in- 
tegral part of the residual shell-model interaction. The; 
shell-model Hamiltonian is usually written in the p-p rep- 
resentation, but it also can be transformed to the p-h rep- 
resentation by means of the Pandya transformation [124]. 
This means that the high- J interaction between pairs can 
translate into the low-J interaction in the p-h channel. 
It is only in the mean-field theory that the division into 
"particle-hole" and "particle-particle" channels appears 
naturally. One way of translating the shell-model results 
into mean-field language is by means of correlators, such 
as the number of T=0 and T=l pairs in the shell-model 
wave function, [69, 112, 120, 125]. 

The extension of the Interacting Boson Models (IBM) 
to the case of pn bosons had to wait until 1980, when 
IBM-3 (only T=l pairs [126]) and IBM-4 [both (T=l, 
S^O) and (T=0, S^l) bosons [127]] were proposed. 
For recent applications of various algebraic models, see 
Refs. [26, 128-138]. 

An alternative strategy to the pn pairing problem is via 
the mean-field approach. Here, the major conceptional 
step was the proposition that quasiparticles are mixtures 
not only particles and holes but also protons and neu- 
trons. The resulting HFB quasi-particle vacuum is a 
superposition of wave functions corresponding to even- 
even and odd-odd nuclei with different particle numbers. 
Unlike in the standard nn and pp pairing cases, the coef- 
ficients of the Bogoliubov transformation are, in general, 
complex. Generalized Bogoliubov transformation, gener- 
alized gap equations, and pn pairing fields are discussed 
in Refs. [3, 53, 55, 57, 62, 65, 67, 68, 101, 139-161]. 

The problem of the spontaneous isospin breaking in 

the mean-field theory was realized soon after the develop- 
ment of the generalized quasiparticle approach [48, 144, 
148] . The symmetry is broken by the independent (sepa- 
rate) treatment of T=l proton and neutron pairing cor- 
relations and by the BCS quasiparticle mean field (the 
generalized product wave function is not an eigenstate 
of isospin). Several techniques have been developed to 
restore isospin. They include the Generator Coordinate 
Method, RPA, Kamlah expansion, iso-cranking, and ex- 
act projection [15, 16, 27, 37, 48, 67, 144, 148, 162-166]. 
It is fair to say, however, that in spite of many attempts 
to extend the quasiparticle approach to incorporate the 
effect of pn correlations, no symmetry-unrestricted mean- 
field calculations of pn pairing, based on realistic effective 
interaction and the isospin- conserving formalism have 
been carried out. 
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III. DENSITY MATRICES IN THE ISOSPIN 

SPACE 

Wc begin with the discussion of the building blocks 
of the HFB theory: one-body density matrices. In the 
HFB theory, expectation values of all observables and, in 
particular, of the nuclear Hamiltonian can be expressed 
as functionals of the density matrix p and the pairing 
tensor k defined as [167] 

p{rst,r's't') = {^\a+,^,t,arst\^), (la) 
k{rst,r's't') = (*|ar.'s't'arst|*), (lb) 

where a^^^ and a^st create and annihilate, respectively, 

nucleons at point r, spin ,s=±^, and isospin t=±^, while 
is the HFB indepcndent-quasiparticle state. Instead 
of using the antisymmetric pairing tensor it is more con- 
venient to introduce the p-p density matrices that can 
be defined in two forms, p or p, denoted by "tilde" and 
"breve" , respectively: 

p{rst,r's't') = -2s'(*|aw-s't'ar.t|*), (2a) 
P{rstys't') = As't'{^ar'-s'-t'arstYi>). (2b) 

In Ref. [5] , p-p density matrix p was used to treat the n-n 
and p-p pairing correlations without the proton-neutron 
mixing. It was then shown that for conserved time- 
reversal symmetry p is hermitian, and leads to p-p lo- 
cal densities that have the structure wliic:li is analogous 
to that of the p-h local densities. However, in the case 
of the proton-neutron mixing studied here, we decided 
to use the p-p density matrix p, because it allows a more 
transparent treatment of the isoscalar and isovector pair- 
ing channels. Detailed discussion of this point will be 
presented in Sec. Ill C below. 

With each of density matrices of Eqs. (la) and (2b) 
three other matrices are associated: 

- the hermitian conjugate matrices defined as: 

p+(rsi,rVi') = p* {r' s't' ,rst), (3a) 
l}+{rst,r's't') = P*{r's't',rst) (3b) 

- the time-reversed matrices defined as: 

p^{rst,r's't') = 4ss'p*{r -st,r' -s't'), (4a) 
p^{rst,r's't') = Ass'p*{r-st,r' -s't'), (4b) 



- the charge-reversed matrices defined as: 

p^{rst,r's't') = 4tt'p{rs-t,r's' -t'), (5a) 
j!F{rst,r's't') = Att' p{rs -t,r' s -t'), (5b) 



where the asterisk stands for the complex conjugation. 

Here and below we present full sets of expressions even 
in those cases when they could, in principle, be replaced 
by verbal descriptions. We do so in order to avoid pos- 
sible confusion at the expense of a slight increase in the 
length of the paper. Wc think that such an approach is 
highly beneficial to the reader, because in many cases 
small but significant differences appear in expressions 
that otherwise could have seemed analogous to one an- 
other. 

The charge-reversal operation C defined in Eq. (5) ex- 
changes the neutron and proton charges, or equivalently, 
flips their isospin projections. Note that the time reversal 
is antilinear while the charge reversal is a linear opera- 
tion, and that they commute with one another. Sym- 
metries of the density matrices can be conveniently ex- 
pressed in terms of just the hermitian conjugation, and 
time and charge reversals. Namely, it follows from defi- 
nitions (la) and (2b) that 

p+ = P, (6a) 
P+ = -r^, (6b) 

where the superscript TC denotes superposition of the 
time (4) and charge (5) reversals. 

For [^P) being an indepcndcnt-riuasiparticlc state the 
density matrices fulfill the following kinematical condi- 
tions 

p*p = /5»/5^^, (7a) 
p = p»p + p»p+, (7b) 

where • stands for integration over spatial coordinates 
and summation over spin and isospin indices, denoted by 
Eda;, e.g.: 



{p» p){risiti,r2S2t2) = {p» p)ixi,X2) ='^dx p{xi,x)p{x,X2) = J d^r'^p{riSiti,rst)p{rst,r2S2t2), (8) 



where we also abbreviated the space-spin-isospin vari- ables by x={rst}. Equations (7) secure the projectivity 
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of the generalized density matrix: 

where 1 := 5{x — x') := 5{r — r')Sss'Stt' and the unitary 
matrix VV, 

transforms the standard generalized density matrix 'R (cf. 
Ref. [167]) to the "breve" representation. 

When the pairing correlations of only like nucleons are 
taken into account, none but the diagonal (off-diagonal) 
matrix elements of density matrix p (p) in isospin in- 
dices are considered. However, in a general case of pair- 
ing correlations between both, like and unlike nucleons, 
the remaining matrix elements become relevant as well. 
Therefore, in the following subsections we specify the 
spin-isospin structure of the p-h and p-p density matrices 
explicitly. 

A. Non-local densities 

The density matrices in the spin and isospin spaces 
can be expressed as linear combinations of the unity and 

I 

p{rst,r's't') 
j){rst,r's't') 

where nt' = (f-jV , , , f f^, ) and &ss' = , ^1^, , 
are the isospin and spin Pauli matrices, respectively, 
which are accompanied by the corresponding unity ma- 
trices, T^^, ~ 5tt' and ct"^, = Sgs' ■ The density matri- 
ces defined in Eqs. (la) and (2b) are now expressed by 
several functions of the pair of position vectors r and 
r' . To avoid confusion, the functions appearing on the 
right-hand sides of Eqs. (11) will be called the (non-local) 
density functions or, simply, densities, unlike the density 
matrices of Eqs. (la) and (2b) appearing on the left-hand 
sides. 



The densities are traces in spin and isospin indices of 
the following combinations of the density and the Pauli 
matrices: 



• scalar densities: 



Pauli matrices. To write the corresponding formulae the 
following notation is assumed. Vectors and vector opera- 
tors in the physical three-dimensional space are denoted 
with boldface symbols, e.g., r or V, and the second rank 
tensors - with sans serif symbols, e.g., J. Scalar products 
of three-dimensional space vectors are, as usual, denoted 
with the central dot: r-V. The components of vectors 
and tensors are labelled with indices a, 6, c and the names 
of axes arc x, y, and z, e.g., r={rx,ry,rz). In order to 
make a clear distinction, vectors in isospace are denoted 
with arrows and scalar products of them — with the cir- 
cle: vow. The components of isovectors are labelled with 
indices i, fc, and the names of iso-axes are 1, 2, and 3, e.g., 
v=(vi,V2,V3). Finally, isoscalars arc marked with sub- 
script "0", and we often combine formulae for isoscalars 
and isovectors by letting the indices run through all the 
four values, e.g., fc=0,l,2,3. 



With this convention the density matrices have the fol- 
lowing form 



(11a) 
(lib) 

I 

- p-h isoscalar and isovector densities: 

Po{r,r') = 5^p(rst,r'si), (12a) 

St 

p{r,r') = Y,Pirst,r'st')i't't, (12b) 

stt' 

- p-p isoscalar and isovector densities: 

po{r,r') = J2^p{rst,r'st), (13a) 

St 

P{r,r') = Y.krstyst'Yn't, (13b) 

stt' 

• vector densities: 

- p-h spin isoscalar and isovector densities: 

so{r,r') = '^p{rst,r's't)&s>s, (14a) 

ss't 

s{r,r') = p{rst,r's't')&s'sh't, (14b) 

ss'tt' 



= jPo{r,r')Sss'6tt' + jSss'p{r,r') o Tt^, + iso(r,r') ■ &ss'Stt' + js{r,r') ■ &ss' o fw , 
= jPo{r,r')das'Stt' + jSss'P{r,r')oftt' + jSo{r,r') ■ &ss'^tt' + js{r,r') ■ &ss' ofw, 
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p-p spin isoscalar and isovector densities: 

so{r,r') = ^p{rst,r's't)&s's, (15a) 

ss't 

g(r,r') = ^ ^{rst,r's't')&s's^t't- (15b) 



ss'tV 



Since the p-h density matrix and the Pauli matrices 
are both hermitian, all the p-h densities are hermitian 
too, 



Po{r,r') 
p{r, r') 

s(r, r') 



p*{r',r), 
s*{r',r), 



(16a) 

(16b) 
(16c) 
(16d) 



and hence, their real parts are symmetric, while the imag- 
inary parts are antisymmetric, with respect to transpo- 
sition of spatial arguments r and r'. 

On the other hand, the unity matrices 0-"^, = 5ss' 
and fj°, = Stt' (scalar and isoscalar) are TC-symmetric 
while the vector and isovector Pauli matrices are TC- 
antisymmetric, i.e.. 



— s —s' 1 



(17a) 
(17b) 



Wc should stress here again that operation TC is an- 
tilinear, and therefore, complex conjugation appears in 
all right-hand-sides of Eqs. (17), although only the Pauli 
matrices Gy and T2 are imaginary. 

Since the p-p density matrix transforms under TC as in 
Eq. (6b), the p-p densities are either symmetric (scalar- 
isovector and vector-isoscalar) or antisymmetric (scalar- 
isoscalar and vector-isovector) under the transposition of 
their arguments, namely: 



po{r,r') = 

p{r,r') = 

so{r,r') = 

s{r,r') = 



-po{r',r), 
p{r',r), 

—s{r', r). 



(18a) 
(18b) 
(18c) 
(18d) 



Equations (16) and (18) are fulfilled independently of any 
other symmetries conserved by the system; they result 
from general properties (6) of density matrices p and p. 



having one spatial argument to distinguish them from 
the non-local densities that have two. Moreover, for lo- 
cal densities the spatial argument will often be omitted 
in order to lighten the notation. 

Following the standard definitions [171, 172], a number 
of local densities are introduced: 

• scalar densities: 



particle and pairing densities: 



Pk{r) 



Pk{r,r')r=r', 
p{r,r')r=r', 



— p-h and p-p kinetic densities: 

nir) = [(V- V')pfe(r,r')],^,„ 
f(r) = [(V.V')p"(r,r')]_,„ 

vector densities: 

- p-h and p-p spin (pscudovector) densities: 



(19a) 
(19b) 



(20a) 
(20b) 



Sk{r) = Sfc(r,r')r 
So(r) = so{r,r') 



r' ) 
r—r' ) 



(21a) 
(21b) 



- p-h and p-p spin-kinetic (pseudovector) densities: 

n{r) = [(V.V')sfe(r,r')]_,„ (22a) 
To(r) = [(V.V')5"o(r,r')]_,„ (22b) 

- p-h and p-p current (vector) densities: 

hir) = ii[(V- V')pfc(r,r')]_,„ (23a) 
Mr) = 5i[C^-'^')/3o(r,r')],^,„ (23b) 

- p-h and p-p tensor- kinetic (pseudovector) densities: 

Fk{r) = i[(V®V'+V'®V)-Sfc(r,r')],^,„(24a) 
Fo{r) = i[(V®V'+V'0V).so(r,r')]_,„(24b) 

• tensor densities: 

- p-h and p-p spin-current (pseudotensor) densities: 



B. Local densities 

In the HFB theory with the zero-range Skyrme inter- 
action [168, 169], or in the local density approximation 
(LDA) (cf. Refs. [167, 170]), the energy functional de- 
pends only on local densities, and on local densities built 
from derivatives up to the second order. These local den- 
sities arc obtained by setting r'=r in Eqs. (12)-(15) after 
the derivatives are performed. They will be denoted by 



Mr) = Jj[(V- V')0^fe(r,r')],^,„ (25a) 
J>) = A[(V-V')®l(r,r')],^,„ (25b) 

where fc=0,l,2,3, and (g) stands for the tensor product of 
vectors in the physical space, e.g., {v'Siw)ab = VaWb and 
[{v<Si w)-z]a = Va{w-z). Note that for particle, pairing, ki- 
netic, spin, spin-kinetic, and tensor-kinetic densities only 
the symmetric non-local densities contribute, while for 
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the current and spin-current densities only antisymmet- 
ric ones contribute. It is then clear that for each p-h den- 
sity there exist both isoscalar and isovector component, 
while for the p-p densities, the isovector component exists 
only for the pairing, kinetic, and spin-current densities, 
while the isoscalar one exists only for spin, spin-kinetic, 
tensor-kinetic, and current densities. 

We note here in passing that the complete list of all 
local densities (up to the derivatives of the second or- 
der) also includes the kinetic and spin-kinetic densities 
in which the two derivatives are coupled to a tensor, i.e., 
VigjV. The resulting local densities arc usually disre- 
garded, because they do not have counterparts to form 
useful terms in the local energy density. There is one set 
of exceptions, which has been overlooked in the system- 
atic construction presented in Ref. [173], and appears in 
the averaging of a zero-range tensor force [171], namely, 
the set of the tensor-kinetic local densities (24). In Sec. 
IV we define terms in the energy density that depend on 
the tensor-kinetic densities. 

All tensor densities (25) can be decomposed into trace, 
antisymmetric, and symmetric parts, giving the standard 
pscudoscalar, vector, and pseudotcnsor components that 
we show here to fix the notation: 

Jk{r) = ^ haair), (26a) 

a=x,y^z 

J(r) = j'-(^)' (26b) 



Jka{r) = ^ eabchbc{r), (27a) 

^a(^) = ^o.bc\c{r), (27b) 



L(r) = \X^{r) + \\^{r)-\j{r)5ab, (28b) 

where /;;=0,1,2,3. 

It follows from Eqs. (16) and (18) that the p-h den- 
sities are all real whereas the p-p densities are in gen- 
eral complex and thus the complex-conjugate densities 
are relevant. The p-p densities become real or imaginary 
only when the time-reversal symmetry is conserved, see 
Sec. VII. 

Instead of the isoscalar and the third component of 
isovector p-h density one can always use the neutron and 
the proton one, e.g., 

Pn{r) = l{po{r)+p,{r)). (29a) 
Pp{r) = i(Po(r)-p3(r)), (29b) 

and just the same for all other p-h densities. Similarly, 
instead of the k=l,2 isovector p-p densities one can use 



the neutron and proton pairing density, i.e., 

Mr) = i(PiW+ip2(r)), (30a) 
Pp{r) = hiPiir)-ip2{r)), (30b) 

and just the same for all other p-p densities. 



C. The TC symmetry 

When constructing the energy density functional (Sec. 
IV) from the local densities (19)-(25) one should ensure 
that it is invariant with respect to: 1° spatial rotations, 
2° isospin rotations, 3° space inversion, and 4° time re- 
versal. All the local densities of Sec. Ill B have definite 
transformation properties with respect to the first three 
of those, l°-3°, so one can easily construct the corre- 
sponding invariants by multiplying densities of the same 
type by one another. For example, a product of any 
pseudovector-isoscalar density with itself, or with any 
other pseudovector-isoscalar density, is an invariant. 

The time-reversal symmetry cannot be immediately 
treated on the same footing, because the time-reversal 
and the isospin rotations do not commute. However, as 
noted in Ref. [174], for problems involving the isospin 
symmetry it is more convenient to use the TC symmetry 
instead of the time-reversal. Indeed, since the charge- 
reversal C is equivalent to a rotations by the angle tt 
about the iso-axis k~2, for conserved isospin the con- 
servation of TC is equivalent to conservation of T alone. 
Therefore, in order to construct the energy density which 
is also time-reversal invariant, we should classify the local 
densities according to the TC symmetry and then mul- 
tiply by one another only densities with the same TC 
transformation properties. 

To this end, we split the p-h and p-p density matrices 
into parts that are symmetric and antisymmetric with 
respect to the TC reversal, i.e., explicitly, 

p±{x,x') = ^(^p{x,x') ±16ss'tt'p*{x,x)j, (31a) 
^±{x,x') = i(/9(x,a:')±16ss'tt'3*(x,x')), (31b) 

where we used a short-hand notation of x={r,~s,—t}. 
In conjunction with the TC transformation properties of 
the Pauli matrices (17), one than immediately obtains 
that the corresponding non-local densities of Sec. Ill A 
are either real or imaginary, i.e.. 



Po±ir,r') 

p±{r,r') 
So±{r,r') 
s±{r,r') 



±Po±iry), 
TpUr,r'), 

±sl{r,r'), 



(32a) 

(32b) 
(32c) 
(32d) 
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and 



Po±{r,r') 
p±{r,r') 
so±{r,r') 

s±{r,r') 



TSo±{r,r'), 



(33a) 

(33b) 
(33c) 

(33d) 



This result shows that real and imaginary parts of non- 
local densities (12)-(15) have opposite TC transforma- 
tion properties. Prom Eqs. (32) one then obtains classifi- 
cation of local p-h densities, namely, the isoscalar densi- 
ties poir), To{r), and Jo(t") are TC symmetric, and So(r), 
To(r), FQ{r), and jo{r) are TC antisymmetric, while the 
isovector densities p{r), T(r), and J(r) arc TC antisym- 
metric, and s(r), T(r), F{r), and j(r) are TC symmet- 
ric. 

The rules of constructing the p-h energy density are 



thus identical to those valid in the case of no proton- 
neutron mixing [172]. On the other hand, from Eqs. (33) 
one obtains classification of local p-p densities, namely, 
real parts of all p-p densities are TC antisymmetric and 
imaginary parts are TC symmetric. The p-p energy den- 
sity must therefore be built by multiplying real parts of 
different densities with one another, and separately imag- 
inary parts also with one another. These rules are at the 
base of the energy density functional constructed in Sec. 
IV. 



IV. THE ENERGY DENSITY FUNCTIONAL 

In the HFB theory the expectation value of Hamilto- 
nian in state |^') is a functional of the density matrices, 
and reads 



EuFB = i^flnm = H[p, p, 3+] = Tr (f . + iTr (f . p + f . 3+) , 



(34) 



where Tr denotes integration over spatial coordinates and summation over spin and isospin indices. Nuclear many- 
body Hamiltonian H, 



H ='^dx'dxf{x',x)a^,ax + i^da;'ida;2da;ida;2 V{x[x'2,xiX2)a^,^a'^,^ax2axi, 



(35) 



r 



is composed of one-body kinetic energy T and two-body 
interaction V, being expressed in (35) by matrix T{x', x) 
and the antisymmetrized matrix elements V{x[x2, X1X2), 

respectively. Matrices F and T are the single-particle 
(p-h) and pairing (p-p) self-consistent potentials, respec- 
tively, 



r{x[,xi) =^da;2da;2 Vp.h(xia;2,a;ia;2)/5(a;2,a;2), (36) 

t{x[,x'2) =^da;ida;2F2Vp.p(a;iS2,a;i52)p(a;i,a;2), (37) 

where ^2=882.52*2^2 and x={r,—s,~t}. In Eqs. (36) and 
(37) we have indicated that the p-h and p-p potentials 
can be determined by different two-body interactions, 
Vp-h and V^-p, called eff'ective interactions in the p-h and 
the p-p channel, respectively. This places further deriva- 
tions in the framework of the energy-density formalism 
that is not based on a definite Hamiltonian (35). More- 
over, effective interactions, Vp_h and Vp.p, are supposed 
to be, in general, density-dependent. 

In the case of the Skyrme effective interaction, as well 
as in the framework of the LDA, the energy functional of 
Eq. (34) is a three-dimensional spatial integral, 



H = / dVH(r), 



(38) 



of local energy density W(r) that is a real, scalar, time- 
even, and isoscalar function of local densities and their 
first and second derivatives. (Isospin-breaking terms, like 
those resulting from different neutron and proton masses 
and from the Coulomb interaction, can be easily added 
and, for simplicity, are not considered in the present 
study.) In the case of no proton-neutron mixing, the 
construction of the most general energy density that is 
quadratic in one-body local densities was presented in 
detail in Ref. [173]. With the proton- neutron mixing in- 
cluded, the construction can be performed analogically 
by including the additional non-zero local densities de- 
rived in Sec. III. Then the energy density can be written 
in the following form: 



nr) = ^ro(r) + ^ (xt(r) + Xt(r)) 

t=o,i 



(39) 



where we assumed that the neutron and proton masses 
are equal. 

The p-h and p-p interaction energy densities, Xt(^) 

and xt, for t=0 depend quadratically on the isoscalar 
densities, and for t=l - on the isovector ones. Based 
on general rules of constructing the energy density, Sec. 
Ill C, one obtains 
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Xo(r) = C^pI + C^'poApo + CJpoTo + C^'^J^ + C^'J^ + C^^:S> + C^'po^ ■ Jo 

+ C^sl + C^so ■ Aso + C^so ■ To + CIjI + C^^ so ■ (V x jo) + C^' (V -80? + C^so ■ -Fo, (40a) 

XI (r) = Cfp^ + c-f V° Ap + ClpoT + Cl'^P + C'l^P + Cff + Cf-^po V ■ J 

+ C^s^ + C^'s- oAs + Cjs- oT+ Cij'^ + CY^s- o (V x j) + CY'{V ■ s)^ + C[s- o F, (40b) 

where x stands for the vector product, and 

Mr) = C^\sof + Ci^'^sS ■ Aso) + C^^s^ ■ To) 

+ Ci\jof + C^m{sS ■ (V X jo)) + C^'\V ■ sof + C^?fi{SS ■ h), (41a) 

XI (r) = C^l^f + C^P^ifi* o + (7[SR(r o ?) 

+ C^Vf + C('\3f + (7f + CY'^r o V ■ 5). (41b) 



In Eqs. (40) and (41) squares always denote total lengths 
in space and/or iso-space, for complex densities taken in 

the complex sense, e.g., |J(r)p = ^ak^lk^ak- In the 

p-p energy density (41) we show only terms in which the 
products of real parts are added to products of imaginary 
parts. According to the rules based on the TC-symmetry, 
Sec. IIIC, similar terms with both products subtracted 
from one another are also allowed. We do not show them 
explicitly, because they have exactly the form of Eq. (41), 
but without complex conjugations and with absolute val- 
ues replaced by real parts of products. 

When the effective interaction is density-dependent all 
coupling constants, C"s and C"s, may also depend on 
density. If this is the case, however, terms that can be 
transformed into one another by integration by parts arc 
not anymore equivalent. Then, five more types of terms 
may appear in the energy density, see Ref. [173]; we do 
not consider such a possibility in the present study. Note 
that in the p-h channel all coupling constants appear in 
two flavors, for t=0 and 1, while in the p-p channel each 
one appears exclusively either for t=0, or for t=l. 

The expression (39) is fairly general. In particular, 
it is not based on any particular two-body interaction. 
However, if one assumes that the underlying two-body 
potential is local and momentum-independent, the form 
of (39) can be simplified and the number of coupling con- 
stants can be reduced. Two particular cases of practical 
interest are discussed in the following. 



A. Local gauge invariance 

Under a local gauge transformation [175], many body 
wave function is multiplied by position-dependent phase 
factor 

A 

|*')=exp{i5^<^(r,)}|*), (42) 



I 

which induces the following gauge transformations of 
density matrices (la) and (2b), 

p'{rst,r's't') = e'^('')-''^(''')/5(rsi,r's'i'), (43a) 
P'{rstys't') = e'^('')+'^('-')3(rst,r'sY). (43b) 

The Galilean transformation is a local gauge transfor- 
mation for (j){r)=p ■ r, where p is a constant boost mo- 
mentum. In analogy to that, one can introduce the local 
momentum field defined by 

p{r) = Vcpir). (44) 

Local and momentum-independent interaction is in- 
variant with respect to local gauge transformation, and 
hence energy densities (40) and (41) must then also be in- 
dependent of the local gauge. The question whether it is 
possible to model nuclear effective interactions in the p-h 
and p-p channels by a local and momentum-independent 
interaction, is open. Therefore, gauge transformation of 
the energy density can, in principle, be respected or not, 
depending on a choice of dynamics one makes. 

It is easy to tell when the local energy densities (40) 
and (41) are local-gauge invariant, because properties of 
local densities (19)-(25) under gauge transformation read 
explicitly 



p'k 


= Pk 




(45a) 


1 


= Tk+2p-jk +P^Pk, 




(45b) 


s'k 


= Sk, 




(45c) 


n 


= Tk + ^p-h+P^Sk- 




(45d) 


j'k 


= jk +PPk, 




(45e) 


F'k 


= Fk+pJk+^k-P + 


p{p-Sk), 


(45f) 


^k 


= h+P<»Sk, 




(45g) 
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where fc=0,l,2,3, and 



where fc=0,l,2,3, and 







(46a) 


f 




(46b) 






(46c) 






(46d) 


% 




(46e) 










+ i(V(8)So)-p-(p-so)p) 


, (46f) 


J' 


= e^'^l 


(46g) 



Since all local p-p densities (46) are multiplied under the 
gauge transformation by phase factors e^"^*^''\ products 
of local p-p densities are not gauge invariant. Therefore, 
all terms not shown explicitly in the p-p energy density 
[see discussion below Eq. (41)] violate the gauge invari- 
ance. On the other hand, products of complex conjugate 
p-p densities and p-p densities may be gauge invariant. 
This obviously is the case for the pairing, spin, current, 
and spin-current p-p densities, while only specific combi- 
nations of kinetic, spin-kinetic, and tensor-kinetic densi- 
ties are gauge invariant. 

Complete list of all p-h and p-p gauge-invariant com- 
binations of local densities reads 



GUr) 
Gl{r) 

Glir) 



= PkTk - Jfc, 

= Sk-Tk- Si 



Sk ■ Tk - 
Pfe V • Jk 



1 72 



1 r2 



ifc> 



Sfe • (V X jk), 

2 



(47a) 

(47b) 
(47c) 



ikba 



Sk ■ Fk - \[^^ikaa^ — ^ JfcabJA 

a ab 

Sk-Fk-lJl + \Jl~\Si, (47d) 



Gl{r) = ■ To) - {^sS ■ Aso) , (48a) 

(r) = SR(so* ■ h) + il V • sol', (48b) 
Gl{r) = 3?(4* f,) - ^p^Apk), (48c) 

where A;=l,2,3. Note that terms of the p-p energy density 

that depend on V x Jq and V • J are not gauge invariant. 

Finally, energy density given by Eqs. (40) and (41) is 
gauge invariant provided the coupling constants fulfill the 
following constraints. 



Ci = 

Of = 

Ci' = 

CP = 

c7' = 



iriT 

T 



for t=0,l, and 



2riF 
\riF 

, 



As 




\riT 

0, 

1 f^r 

0. 



(49a) 
(49b) 
(49c) 
(49d) 
(49e) 

(50a) 
(50b) 
(50c) 
(50d) 
(50e) 



B. Skyrme interaction energy functional 

The Skyrme interaction [168, 169] is a zero-range lo- 
cal force that depends on relative momenta up to the 
second-order. The complete list of terms giving its ma- 
trix element in the position-spin-isospin representation, 
including the tensor components [171, 176], reads 



y(r;s'ii'ir^s'24,risiii v^s^X^) = [t^i}" ^ x^P") + It^" + xzP^)p'^ (^(ri +r2)) 



+ iii(5^+a;iP<^) 



-r 2''e 



fe'* • S • fe'* -h fc • S • fc 



t2(^'' + X2p'')k'* ■ k + tok'* ■Sk + iWoS • 



k'* X A; 



pa pT pM 



(51) 

)^12, 



r 



where 



and 

p-r 



<5fifi<5f^f2, 



(52a) 
(52b) 



are the spin and isospin unity and exchange operators, 
respectively, and 



\'yK,s'2s,s2+^s',s^ •o-,^s2)= '5,is2'54.i,(53a) 
liK'^t'^t^t^ + ^ti*i ° ^t^ts) = <5fif2<5t^fi> (53b) 



— ^(rr°- 

2 I'^s'iSi 



±.b 

S'S2 



s'iSl"s2S2/ 



(54a) 



(54b) 
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are two-body vector and tensor spin operators, respec- 
tively. The relative momentum operators, 

^ = 5i(Vi-V2), (55a) 
k' = 5l(Vi-V^), (55b) 

act on the delta functions in 612, 

Si2{r[r2,rir2) = S{r[-ri)5{r'2-r2)S{ri-r2) 

= S{r[—r2)6{r'2—ri)5{r2 — ri). (56) 

This action has to be understood in the standard sense 
of derivatives of distributions. 

Whenever the Skyrme interaction (51) is inserted into 
integrals, like in Eqs. (35)-(37), the integration by parts 
transfers the derivatives onto appropriate variables in the 
remaining parts of integrands. 

Numbers are equal to +1 or —1 depending on 
whether in a given term the power of momentum k is 
even or odd, respectively. Skyrme interaction written in 
the form of the integral kernel (51) is explicitly antisym- 
metric with respect to exchanging left or right pairs of 
variables pertaining to particles 1 and 2. 

The Skyrme HFB energy density can be calculated 
by inserting the Skyrme interaction (51) directly into 
expressions (36), (37), and (34). Results for the p-h 
channel were published by many authors, see, e.g., Refs. 
[169, 172, 175, 177], although often some terms of inter- 
action (51) were neglected and/or restricted symmetries 
were used. Results for the p-p channel were previously 
published with tensor terms and the proton- neutron mix- 
ing neglected [5] . Here we aim at presenting the complete 
set of results. 



J 



Calculations leading to expressions for the Skyrme en- 
ergy density are tedious, but can be efficiently performed 
by noting two simplifying facts. First, the two-body spin 
operators obey conditions. 



SP" = S, 
SP" = S, 



(57) 
(58) 



and hence only terms up to linear in spin and isospin 
Pauli matrices appear in the antisymmetrized interac- 
tion. Second, the Pauli matrices in (51) pertain to the 
p-h coupling channel, while the momenta - to the p-p 
coupling channel. Hence, calculations may become very 
easy once a common, p-h or p-p, coiipling channel is used 
for all the elements of interaction. This requires either 
recoupling momenta to the p-h channel, or recoupling 
the Pauli matrices to the p-p channel. To this end, we 
separately consider the p-h end p-p energy densities. 



1. The p-h channel 

In the p-h energy density, indices of the Pauli matrices 

are contracted directly with density matrices of particles 
1 and 2, and immediately give non-local densities through 
appropriate traces in Eqs. (12) (15). However, the rel- 
ative momentum operators (55) affect both particles at 
the same time, and hence have to be first recoupled to 
forms where the two particles are acted upon indepen- 
dently, i.e.. 



k'* k 
k'* X k 
k'* ®k'* +k(S:k 



k'* (g)k + k(g) k' 



l{Kf + Kl -Ki K2- 4fci • fca) + 3(Vi • V[ + V2 • V^), 

X (fca - kl), 

\{ki®k2 + k2®ki) -\{kx 
\{v[ ® Vi + V2 ® V2), 



K2) 

1 •i9 Kl 

\{^i®V'i 

\{ki®k2 



K2 ® K2) 
V2 ® v'2) 
k2 kl) 
V2 ® V2) 



\{ki (^k2 + k2® kl) 
i(V'i0Vi-h V20V2) 



^2+^2® kl) 



(59a) 
(59b) 
(59c) 

(59d) 

(59e) 



where 



and 



kl 
k 

kl 
k2 



ii(Vi-v;) 

A(V2-V^). 



-*(Vi 
-i (V2 



V'l), 



(60a) 
(60b) 

(61a) 
(61b) 



Final results can now be easily obtained by noting that 
relative momenta (60) lead to the current densities (23a) 
and (25a), total momenta (61) lead to derivatives of local 
densities, and the scalar and tensor products of individual 
momenta lead to kinetic densities (20a), (22a), and (24a). 

The zero-order (density-dependent) p-h coupling con- 
stants of the energy density (40) are expressed by the 
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TABLE I: Second-order coupling constants of the p-h energy 
density (40) as functions of parameters of the Skyrme inter- 
action (51), expressed by the formula: C = ^{ati + btixi + 
ct2 + dt2X2 + ete + fto + qWo). 
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Skyrme force parameters as 

Co' = 1*0 + ii3P?(r), (62a) 

Q = iio(2aro - 1) + ^ts{2x3 - l)p?(r), (62b) 

Cf = -ito(2xo + l)-^i3(2a;3 + l)p?(r), (62c) 

Ct = -5*0 - 35*3P^W. (62d) 

and the second-order coupling constants are given in Ta- 
ble I. One can immediately see that the gauge-invariance 



conditions (49) arc fulfilled. This is so because the 
momentum-dependent terms of the Skyrme interaction 
obey the Galilean invariance [172, 175] 

Since seven Skyrme force parameters define twenty 
four second-order p-h coupling constants, in the result- 
ing Skyrme energy density there is a high degree of de- 
pendency. First, as is well-known [178], a single spin- 
orbit parameter Wq determines all four spin-orbit cou- 
pling constants C^'^ and C^'' , for t=0 and 1. Second, 
four Skyrme parameters, ti, .xi, 1,2, and X2, uniquely de- 
termine four coupling constants Cf^'' and C^, for t=0 
and 1. Third, two tensor Skyrme parameters, te and to, 
uniquely determine either isoscalar or isovector coupling 
constants, C^'^ and Cf. Once such a role of the seven 
Skyrme parameters is fixed, values of the remaining cou- 
pling constants are also uniquely fixed. 

2. The p-p channel 

In the p-p energy density, each operator of the relative 
momentum, k' and k, acts on variables of the same den- 
sity matrix, and thus no recoupling is necessary. Terms 
of the interaction that are linear in momenta then lead to 
current densities (23b) and (25b), while terms that are 
quadratic in momenta lead to derivatives of local den- 
sities and to kinetic densities (20b), (22b), and (24b), 
because 

fe' = -i(Vi + V2)' + Vi- V2, (63a) 
fe0fe = -^(Vi -h V2) O (Vi -h V2) 

+|(Vi V2 + V2O Vi). (63b) 

However, in the p-p energy density, indices of Pauli 
matrices couple together the two density matrices, and 
hence do require recoupling to the p-p channel. These re- 
coupling formulae can be obtained by means of the stan- 
dard algebra of angular momentum. A sum of the three 
Clebsch-Gordan coefficients appropriate to the present 
case reads [179] 



4S2S2 ^ (5SiAi/xi|isi)(5 -S2A2/X2I5 -s'2){Xi -/X1A2 -M2|AAt) 
^ll^t2 

= ^ (-l)^^-^=+''(2r + l)(2Z + l) I I I xA X {^s'J'm'\^s[){^silm\^S2){l' ~rn'l-m\Xti) (64) 

l'm',lm [ Z' Z A J 

Taking relevant combinations of Ai, A2 = 0, 1 and A = 0, 1, 2, one obtains: 

44s2^^;,_4_,,,._,,, = ^S,>^s[Ss,s,+^(r;,^s[-'^s2s,, (65a) 

44«2A",,-4,,i,-s. = -i<54.'/,,,, -(T,,,,, (65b) 

4s2'S2^si,-4si,-s2 = X , (65c) 

4452S:? = -f + <,J + Saba;,,, ■ a,,,, (65d) 
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TABLE 11: Second-order coupling constants of the p-p energy 

density (41) as functions of parameters of the Skyrme inter- 
action (51), expressed by the formula: C = ^{ati + ht\Xi + 
Ct2 + dt'2X2 + ete + fto + gWo)- 
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and the formulae similar to Eqs. (65a) and (65b) are ob- 
tained for 5'^ and P^, respectively. 

The two zero-order (density-dependent) p-p coupling 

constants of the energy density (41) are related to the 
Skyrme parameters in the following way: 

C-o^ - lto{l + xo) + i^h{l + x^)p'^{r), (66a) 
CI = lto{l-xo) + ^h{l-Xi)p'^(r), (66b) 

and the second-order p-p coupling constants arc given in 
Table 11. Similar to the p-h case, the gauge-invariance 
conditions (50) are met. 

Equivalcntly, the density-dependent zero-range pairing 
force Vpair can be used in the p-p channel [180-183], 



Vpair(r,r') = /pair(r)5(r - r') 



(67) 



for 



/pair(r) = Vo\l+ XoP" 



Po{r) 



(1 + XsP' 



(68) 

where P°' is the spin-exchange operator (53a). In such a 
case, coupling constants (66) read 



= iyo(l + a;o)-|Vo(l + a;3) 
= iyo(l-xo)-iyo(l-X3) 



Po{r) 



Pc 

Po(r) 



(69a) 
(69b) 



Note that when only the isovector pairing is used, as in 
most LDA applications to date, the exchange parameters 
Xq and 2:3 are redundant in the definition of the isovec- 
tor coupling constant , and hence are usually set to 
0. However, if one wants to independently model the 
isoscalar and isovector pairing intensity, one has to use 
non-zero values of xo and xs. 



For the Gogny interaction [167], the zero-range densi- 
ty-dependent term with a=l/3 was used in order to 
enforce proper saturation properties. The correspond- 
ing exchange parameter 0:3=1 was used to prevent this 
zero-range force from contributing to the isovector pair- 
ing channel. However, such a choice, when applied liter- 
ally to the proton-neutron mixing case, might lead to a 
very strong repulsive isoscalar pairing interaction. 

The term of x coming from the spin-orbit interaction 
contains the combination of components of the p-p spin- 
current density J, 



ab 



bb 



Jab ° Jfea 



(70) 



that is different from that coming from the tensor to term, 



ab 



■i*ab ° Jab - |Joa ° Jbb " fJab ° Jba) 

= _||/|2+5|5|2_1|5|2^ (71) 



and from that coming from the central t2 term. 



ab 



Therefore, by setting appropriate values of the t2{l + X2), 
Wq, and to parameters, one can obtain arbitrary values 
of the spin-current coupling constants Cf, Cf, and 
cf. Similarly, parameter t2{l — X2) allows for fixing 
an arbitrary value of the current coupling constant Cq. 
On the other hand, parameter ti{l + xi) defines two 
isoscalar coupling constants, (7,^* and Cj, parameter 
te defines another two isoscalar coupling constants, Cg^* 
and Cq , and parameter ti{l — xi) defines two isovec- 
tor coupling constants, ct' and Cl; hence, these pairs 
of coupling constants are not independent from one an- 
other. These three pairs of dependencies reflect, in fact, 
the three gauge invariance conditions (48). In this way, 
seven Skyrme force parameters determine ten coupling 
constants in the p-p channel. Finally, the Skyrme inter- 
action does not give any non-zero values for the spin-orbit 
coupling constants Cq-' and Cf'' ■ Therefore, up to the 
gauge invariance conditions, the Skyrme interaction fully 
determines the energy density in the p-p channel. 



V. THE P-H AND P-P MEAN FIELDS 

By varying the energy functional (34) with respect to 
the density matrices one obtains the p-h and p-p mean 
field Hamiltonians, 
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h{r's't',rst) 



h{r's't',rst) = 



S^+{rst,r's't') 



= t{r's't', rst) + tr{r's't', rst). 



(73b) 



r 



The rearrangement potentials f r and f r result from the 

density dependence of effective interactions on the p-h 
and p-p densities, respectively. Usually effective inter- 
actions are assumed to depend only on the p-h density 
matrix (most often, only on the isoscalar particle density 
po)- In that case the p-p rearrangement potential van- 
ishes. However, one cannot forget that the dependence 
of the p-p interaction on the particle density results in 
a corresponding contribution to the p-h rearrangement 
potential. In what follows, to simplify the presentation 
we do not show the rearrangement terms explicitly. 

Within the LDA, the mean-field Hamiltonians being 
originally, like the Skyrme interaction of Eq. (51), either 
distributions or derivatives of distributions, can, when 
acting as the integral kernels, be expressed as local, mo- 
mentum dependent operators, i.e., 

h{r's't',rst) = 6{r - r')h{r; s't' , st), (74a) 
h{r's't',rst) = S{r - r')h{r;s't',st). (74b) 



The kinetic energy term in Eq. (73a) is already expressed 
in such a form. The mean-fields Hamiltonians are the 
second-order operators in momentum and matrices in the 
spin and isospin spaces. The isospin structure of the local 
p-h and p-p mean-field Hamiltonians reads 



h{r;s't' ,st) = ho{r;s' ,s)5t>t + h{r;s' ,s) o fft, (75a) 
h{r; s't' , st) = ho{r;s',s)dt't + h{r;s',s)oTt,t, (75b) 



respectively. The isoscalar and isovector parts of the p-h 
mean-field Hamiltonian can be presented in the compact 
form 



hk{r; s', s) = - — V^Sg'sSko + UkSs>s + • as's + ^ [ik^s's + (Bfe • &s's)] ■ V + ■ [hSs's + (Bfe • &s 



-V • [MkSs's + Ck ■ as's] V - V • Dfe as's ■ V 



(76) 



for fc = 0, 1, 2, 3, and where 

(B ■a)a = Y, Bab^', (77) 

6 

for a = x,y, z, is the a'th component of a space vector. 
The names of symbols arc inspired by those introduced in 
Ref. [172]. Since the p-h density matrix is hermitian, the 
p-h mean-field Hamiltonian is also hermitian and, thus, 
all the potentials, Mk, Uk, B/j, Ck, Dk, Ik, and Sfc are 
real. 

The general form of the mean-field Hamiltonian (76) 

can be constructed from the momentum — i V and spin & 
operators, based only on the symmetry properties. Apart 
from the one-body kinetic energy [the first term in Eq. 
(76)], the expansion in momentum gives: (i) zero-order 
terms with scalar (Uk) and pseudovector (Sfe) potentials, 
(ii) first-order terms with vector (Ik) and pscudoten- 
sor (B;j) potentials, (iii) second-order-scalar terms with 



scalar (Mk) and pscudoscalar (Cfe) effective masses, and 
(iv) second-order-tensor terms. In principle, the most 
general form of the last category should involve tensor 
and third-order-pseudotensor potentials. However, in Eq. 
(76) we show only the particular form of it that corre- 
sponds to the energy density (40). 

According to Eqs. (73) the p-h mean-field Hamiltonian 
is the functional derivative of the energy functional over 
the hermitian p-h density matrix. Functional derivatives 
of integrals of type: 

(7p>= y (^(n -r2)/(ri)p(ri,r2)dVidV2, (78) 

where function / is treated as independent of densities 
and p represents a p-h non-local density, can easily be 
calculated using Eqs. (12) and (14). Bearing in mind 
that 
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5p{riSiti,r2S2t2) ^, i\x a x x 

Sp{rst,r's't') = ^ir,-r)S{r2-r)S,,,S,,,5,,,5,,„ 

(79) 



one has 



'^'^^ ' '>(r'-r)J(r)a...ft,, (80b) 



Sp{rst, r's't') 



for A; = 0,1,2,3. The functional derivatives of integrals of local differential densities are obtained from Eqs. (80) 

through integration by parts. Then, the functional derivatives become dependent on derivatives of the Dirac delta 
function and thus, in accordance with Eqs. (74), again act as local differential operators. They read: 



-^(/Jfe) _ ^^s(r'-r){Vf{r) + f{r)V)5,,,T^,„ (81a) 



Sp{rst, r's't') 



'^^^^'^^ = TAr' - r){VJ{r) + /(r)V„),T^,f,'5„ (81b) 



Sp{rst, r's't') 



Sp{rst,r's't') ~ ^^"^ ^/af[r)Vhds'sT,,t 



S{r'-r) Vafir)Vbds'sTt%, (82a) 



for A; = 0, 1, 2, 3 and a,b,c = x, y, z. Calculations of the functional derivatives over the density matrix are equivalent 
to the rules for variations over single-particle wavefunctions given by Engel et al. [172]. Using formulae given above, 
Eqs. (80)-(82), one obtains the following relations between the potentials defining the p-h mean field (76) and the 
local p-h densities defining the energy density (40), 

(83a) 
(83b) 
(83c) 
(83d) 
(83e) 
(83f) 
(83g) 

for k = 0,1,2,3. All coupling constants Ct in Eqs. (83) are taken with t=Q for /c = (isoscalars), and with t=l 
for k =1,2,3 (isovectors). Symbol 5 is the unit space tensor, and e • J stands for the antisymmetric space tensor 
with components: (e • J)ab = X^c ^acbJc, so that, according to Eq. (77), its action on a vector is obviously the vector 
product: {e ■ J) ■ & — J x &. 

The p-p mean-field Hamiltonian has the following isoscalar and isovector components: 

hoir; s', s) = So • + ^{ V • loS^'s + hSs's • v| - V • [Co • V -V -Do as's ■ V, (84a) 
Mr; s', s) = &{r)6s's + ^{ V • [i(r) • <t,.,] + [i(r) • <t«.J • v} - V • 7^5y«V. (84b) 



Uk{r) 


= 2Cf pfe + 2Cf ''Apfe + ClTk + V • Jk 






Sfe(r) 


= 2C*Sfe + 2(Cf * - Cf')Ask - 2C7"V X 


(V X 


Sk) + CjTk + Cf Ffe + Cj'V X jk 


h{r) 


= 2Cljk + c7'V xsk, 






Bfc(r) 


= 2C/" JfcJ - 2Ci^e ■ Jk + 2C/24 + cJJe 






Mfe(r) 








Ck{r) 


= Sk, 






Dk{r) 


= Cfsk, 







V 



Contrary to the p-h Hamihonian (76), the p-p Hamilto- j)^^ j^^ m, tj, and B are, in general, complex quan- 
nian (84) can be non-hermitian, because potentials Co, 
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titles. This is so, because the p-p density matrix is, in 
general, not hermltian. Therefore, the energy functional 
should be treated as a functional of both p and p+. 

The p-p mean-field Hamiltonian is the functional 
derivative of the energy functional over /5+ , whereas the 
hermitian conjugate Hamiltonian is the functional deriva- 
tive over p. The p-p densities are, according to Eqs. (13) 
and (15), functions of p, while the complex conjugate 



J 



densities are functions of p"*". 

When calculating the p-p functional derivatives, one 
cannot forget that the p-p density matrix fulfills symme- 
try condition (6b), implying that the p-p densities arc 
either symmetric or antisymmetric functions, Eqs. (18). 
Therefore, the calculation of functional derivatives over 
cither p or p^ is similar to that leading to Eqs. (80)-(82), 
however, instead of Eq. (79) one has: 



(5/3+(ri.siti,r2S2^2) 
Sl}+{rst,r's't') 



= 5{ri - r)S{r2 - r')6sisSs2s'StitSt2t' - IQss'tt'Sin - r')d{r2 - r)5si-s'Ss2-sSti-t'St2-t (85) 



In the expressions for functional derivatives, this gives either cancellation or addition of terms coming from the two 
components of the right-hand side of Eq. (85). Finally, the non- vanishing derivatives are 



S{fp*) 



5p+(rst,r's't') 



Sp+{rst,r's't') 



2S{r' - r)f{r)5s'sh't, 
26{r' - r)/(r)(T,.,f°„ 



(86a) 
(86b) 



HfJo) 



6p+{rst,r's't') 



6p+{rst,r's't') 



7-T = -iS{r'-r){Vaf{r)+f{r)Va)K,Tt>t, 

I O'T' \ 



(87a) 
(87b) 



5{fVaV',p*) 

5p+{rst,r's't') 

5j}+{rst,r's't') 



-26{r' -r)Vafir)Vb6,,A't, 



= -25(r'-r)V„/(r)V6<T,^,f°„ 



(88a) 
(88b) 



for a,b,c = x, y, z. 

Using Eqs. (86)-(87) one obtains the following relations between the potentials defining the p-p mean-field Hamil- 
tonian (84) and the local p-p densities defining the energy density (41): 



So(r 
/o(r 
C'o(r 
£>o(r 

tf{r 
i(r 



= + 2{Ct - Cj')Aso - 2Cj'V x (V x sq) 

= 2Cfjo + Cj'V xso, 

= CqSq, 

— (^0 *o, 

= 2CPj+2Ct''l^p + Cl^+CY^V -Jk, 

= 2Cfh - 2C(^e ■ 5 + 2Cpl + CY-^e ■ Vp, 

= cip. 

I 



(89a) 
(89b) 
(89c) 
(89d) 

(89e) 
(89f) 
(89g) 



In the case of the zero-range pairing force (67), the isovec- 
tor p-p potential is proportional to the p-p isovector den- 



sity while the isoscalar field has a very different structure, 
i.e., it is proportional to the scalar product of spin & and 
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the p-p spin density Sq. This immediately suggests that 
there exists a connection between the isoscalar pairing 
and the p-p spin saturation, which is influenced by the 
spin-orbit splitting. In this context, let us remind the 
shell-model study [20] which discusses the relation be- 
tween the magnitude of the T=0 pairing and the spin- 
orbit splitting. 



VI. THE HFB EQUATIONS 

Minimization of the energy functional of Eq. (34) with 
respect to the p-h and p-p density matrices, which fulfill 
Eqs. (7) under auxiliary conditions 



y dVp„(r) = N, 
J dVpp(r) = Z, 
leads to the HFB equation of the form: 



(90a) 
(90b) 



= 0. 



The generalized density matrix TZ is given by Eq. (9) and 
the generalized mean-field Hamiltonian is defined as 



a system of eight second-order differential equations, in 
general with complex coefficients. Usual four dimensions 
corresponding to upper and lower HFB components and 
to two spin projections are here multiplied by another 
factor of two due to the isospin projections. Altogether, 
Eq. (94) corresponds to a system of sixteen equations 
within the domain of real numbers. When specific sym- 
metry conditions are imposed on solutions, this number 
can be reduced in a standard way, sec Ref. [184] for the 
analysis pertaining to spherical symmetry 

The energy spectrum of generalized mean-field Hamil- 
tonian has been discussed in Ref. [5] . The only difference 
with the present case is that here the eigenvalue problems 
for neutrons and protons in Eq. (94) cannot be separated. 

It is well known, that the eigenvalues of H appear in pairs 
of opposite signs. For each quasihole state of energy E 



$(rsi; E) 



( 'first; E) 
\ il){rst;E) 



(96) 



^9]^) there exists a quasiparticle state 



^{rst; -E) = 4st 



i)*{r-s-t; E) 
(p*{r-s-t;E) 



(97) 



H = wnw+ = 



h-X 



(92) 



h+ -h^^ + Xj 

with the Lagrange multiplier given by 

A = i(A„-^Ap)-^i(A„-Ap)f^ (93) 

where A„ and Ap are the neutron and proton Fermi levels, 
respectively. 

The usual method of solving the HFB equation (91) is 
to solve in a self-consistent way the eigenvalue problem. 



n{x',x)»<l>{x;E) = E<i>{x';E), 



(94) 



for the generalized mean-field Hamiltonian, and then to 
construct the generalized density matrix. 



^(a;,x') = ^$(a;;£;)$+(x';i;), 



(95) 



as a projection operator onto the set of the quasihole 
(occupied) states $ belonging to a subset of energy spec- 
trum, £. For a local mean-field Hamiltonian, Eq. (94) is 



belonging to energy —E. In the case of absence of exter- 
nal fields, bound states (when ip and ip are both localized) 
exist only when both Fermi levels, A„ and Ap, are nega- 
tive. Discrete quasihole energy levels lie within the range 
C < E < — £, where C = max(A„, Ap) < 0. The ground- 
state solution corresponds to occupying states having 
negative energies; then the set £ consists of a number 
of discrete levels lying inside segment (£, 0) and the con- 
tinuous spectrum with —oo<E<£. 

Traditionally, one solves Eq. (94) for the quasiparticle 
states of positive energies rather than for the negative 
ones. Then, the discrete spectrum is within the segment 
< E < —C and energies E > —£ belong to the con- 
tinuum. Having found the wavefunctions ^{rst;E) for 
E > one uses Eq. (97) to construct the density matrix. 



I.e., 



'k{x,x') = Y^<^{x;-E)^+{x';-E). (98) 



E>0 



The p-h and p-p density matrices are then expressed as 
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p{rst,r's't') = 16ss'tt'^tp*{r-s-t;E)'ip{r' -s' -t']E), 

E>0 

p{rst,r's't') = 16ss'U'^^*{r-s-t;E)ip{r' -s' -t';E). 

I 



E>0 



(99a) 
(99b) 



VII. CONSERVED SYMMETRIES 

Conserved and broken symmetries are one of the most 
important elements of description of many-body systems. 
Within the mean-field approach, the theorem about self- 
consistent symmetries [167] tells us that mean- field states 
may or may not have all the symmetries of the Hamilto- 
nian, depending on interactions and the system studied. 
Within the HFB approach, the symmetry is conserved 
when the generalized density matrix TZ and the general- 
ized Hamiltonian H both commute with the symmetry 
operator ZY, i.e., ['R.,L{]=0 and ['H,h{]=0, or 



■R, 

n, 



where 



U 

u* 



(100a) 
(100b) 



(101) 



and ?7 is a unitary matrix of the single-particle symme- 
try operator. For the "breve" representation used in the 
present study, the symmetry operator is given) by [cf. 
Eq. (9)] 



U = WUW+ = 



U 

u^^ 



and then 



i{'RLi+ = a, 



(102) 



(103a) 
(103b) 



In the previous sections we have presented the most gen- 
eral set of expressions pertaining to the situation when 
no symmetries were a priori conserved. Below we discuss 
consequences of conserved symmetries. 



A. Proton-neutron symmetry 

The standard case of no proton-neutron mixing can 
be described by the conserved proton-neutron symmetry 
given by 

Upn = i exp(-i7rT3) = i exp(-|7rf3) = fs. (104) 

In other words, the iso-3 signature (multiplied by i) is 
then the conserved symmetry. Note that conservation 



of projection of the isospin on the third axis (the charge 
conservation) would require that the iso-3 rotation about 
an arbitrary angle be conserved, while the iso-3 signature 
corresponds only to rotation about tt. Within the HFB 
approach, the charge symmetry is broken in the same 
way as is the particle number symmetry. 

Since the TC-transformed symmetry operator reads 
fs, we obtain from Eq. (103a) that 



jjTC 
^ pn 



rsPh = -p, 



(105a) 
(105b) 



and analogous properties hold for the mean-field Hamil- 

tonians, h and h, respectively. It is then clear that with- 
out the proton-neutron mixing the p-h density matrices 
and Hamiltonians have only the k=0 and 3 isospin com- 
ponents, while the p-p ones have (in the "breve" repre- 
sentation) only the k=l and 2 isospin components, cf. 
Eqs. (29) and (30). 



B. Time-reversal symmetry 

In the case of time-reversal invariance, p^=p and 
p^=p, see Eqs. (4), the p-h and p-p densities fulfill addi- 
tional conditions, 

po(r,r') = pUr,r'), (106a) 

Pk{r,r') = -(-l)V^(r,r'), (106b) 

So{r,r') = -s*ir,r'), (106c) 

Sk{r,r') = (-l)'=4(r,r') (106d) 

and 

Po(r,r') = pUr,r'), (107a) 

Pk{r,r') = -{-lfpl{r,r'), (107b) 

so(r,r') = -sS(r,r'), (107c) 

h{r,r') = {-!)>' sUry), (107d) 

where A;=l,2,3. Due to the fact that the fc=2 Pauli matrix 
is imaginary, the time reversal acts differently on the 
k=2 isovector components than on the k=l,3 components 
of all isovector densities. At the first sight, this seems 
to be a bizarre property. Indeed, the isospin quantum 
number is introduced to take into account the fact that 
there are two kinds of nucleons in nature, and each kind 
has its own, apparently unrelated to one another, time- 
reversal operation. 
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However, the use of the standard isospin formaUsm im- 
pUes something more. Namely, the neutron wave fune- 
tion (isospin up) can be obtained from the proton wave 
function (isospin down) by an action of the (real) Pauli 
matrix. Therefore, the relative phases of the neutron and 
proton wave functions arc fixed by the phase convention 
that has been used to choose the isospin Pauli matrices. 
As a consequence, the time-reversal properties of neu- 
trons and protons are not any more independent from 
one another. Of course, this is not a spurious quirk of 
the mathematics we use, but a reflection of a deeper fact 
that by mixing the neutron and proton wave functions 
we introduce complex mixing coefficients that do affect 
the time-reversal properties of the mixed wave function. 
Conservation of the time reversal means that these mix- 
ing coefficients must follow rules dictated by the time 
reversal, which implies differences between the k=2 and 
fc=l,3 iso-directions. Therefore, we see here that from 
basic arguments it follows that conservation of the time 
reversal must imply the isospin symmetry breaking. The 
only iso-rotations that are compatible with the time re- 
versal are those about the fc=2 iso-axis. (The influence 
of the time-odd fields on the magnitude of the Wigner 
energy was pointed out in Rcf. [12].) 

Table III summarizes properties of p-h and p-p densi- 
ties under the exchange of their spatial arguments. When 
no conserved symmetry is imposed, all densities are com- 
plex, and their real and imaginary parts are either sym- 
metric or antisymmetric. For conserved time reversal, 
all densities become either real or imaginary, and are 
either symmetric or antisymmetric. Recall that sym- 
metric parts contribute only to particle, kinetic, spin, 
spin-kinetic, and tensor-kinetic local densities, while the 
antisymmetric parts contribute only to the current and 
spin-current local densities. Therefore, local densities 
are complex, real, imaginary, or vanishing, depending on 
whether time-reversal, proton-neutron, or both symme- 
tries are conserved. Table IV presents these properties 
for all local p-h and p-p densities. 

In previous studies, e.g., in Refs. [3, 62, 65, 152, 153], 
the r=l pairing fields were associated with the real part 
of the pairing tensor, while the T=0 pairing was repre- 
sented by the imaginary part of the pairing tensor. Such 
a structure was obtained for specific phase conventions 
and symmetries. On the other hand, as shown in Table 
IV, the general case corresponding to no conserved sym- 
metries (e.g., for rotating states) requires that all the pn 
densities be complex. 

To summarize this subsection, we now enumerate all 
non-zero densities when the time reversal is conserved or 
not, and/or the proton- neutron mixing is present or not. 
By counting as one density each component of a vector, 
tensor, or isovector, we obtain the following four options: 

1° Time reversal broken plus proton-neutron mixing: 

- 23 real p-h isoscalar densities: poi'^'), ''o(r), Jo(''), 
So(r), To(r), jo(r), and Fo(r), 

- 69 real p-h isovector densities: p(r), f{r), J(r), 



TABLE III: Symmetries of the p-h (left) and p-p (right) densi- 
ties in general case (no conserved symmetries imposed), and 
in case of the time-reversal symmetry conserved. Real (5ft) 
and imaginary (9) parts are symmetric (S) or antisymmetric 
(A) under exchange of their spatial arguments, as indicated 
in the Table. 



density 


general 


time-rev. 


density 


general 


time- rev. 


Po(r-,r') 


S 


A 


S 





po(r,r') 


A 


A 


A 


P2(r,r') 


s 


A 





A 


P2(r,r') 


S 


S 


S 


pi,3{r,r') 


s 


A 


s 





pi,-i{r,r') 


S 


S 


S 


so{r, r') 


s 


A 





A 


S()(r,r') 


S 


S 


S 


S2{r, r') 


s 


A 


s 





S2{r,r') 


A 


A 


A 


si,3(r,r') 


s 


A 





A 


si,3(T',r') 


A 


A 


A 



TABLE IV: Properties of the local p-h and p-p densities in 
general case (no conserved symmetries imposed), and in case 
of the time-reversal, proton-neutron, or both symmetries con- 
served. The fc=0,l,2, or 3 isospin components of densities are 
complex (C). real (R). iniagiuary (I), or zero (0), as indicated 
in the Tal)k>. 





general 


time- rev. 


prot. 


-neut. 


both 


k 





1 


2 


3 





1 


2 


3 





1 


2 


3 





1 


2 


3 


pk 


R 


R 
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s{r), f{r), j(r), and F{r), 

- 12 complex p-p isoscalar densities: So{r), Tb(r), 
i'o(r), and Fo(r), 

- 33 complex p-p isovector densities: p{r), f{r), and 
J>), 

2° Time reversal conserved plus proton-neutron mixing: 

11 real p-h isoscalar densities: poir), To{r), and 
Mr), 

- 30 real p-h isovector densities: pi,s{r), Ti^3{r), 
Ji,3(r), S2(r), T2(r), j2{r), and h{r), 

- 12 imaginary p-p isoscalar densities: So{r), Tb(r), 
jo(r), and Fo{r), 
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- 33 p-p isovector densities, 22 real: pi^s{r), Ti_3(r), 
Ji,3(r), and 11 imaginary: P2{r), T2(r), J2(t'), 

3° Time reversal broken, no proton-neutron mixing: 

- 23 real p-h isoscalar densities: po{r), To(r), ^^{r), 
So(r), To(r), jo(r), and Fo(r), 

- 23 real p-h isovector densities: P3{r), T3(r), J3(r), 
S3(r), T3(r), Mr), and F3(r), 

- 22 complex p-p isovector densities: pi,2{r), Ti^2{r), 
and Ji,2(r), 

4° Time reversal conserved, no proton-neutron mixing: 

- 11 real p-h isoscalar densities: Poif), To{r), and 
Mr), 

- 11 real p-h isovector densities: P3{r), T3(r), and 
h{r), 

- 22 p-p isovector densities, 11 real pi{r), fi(r), 
Ji(r), and 11 imaginary M^), M^), "M^)- 

VIII. CONCLUSIONS 

Experimental studies of the heavy N^Z nuclei have 
sparked renewed interest in physics of pn correlations, 
especially pn pairing. While the appearance of the T=l 
pn pairing is a simple consequence of the charge invari- 
ance, in spite of vigorous research, no hard evidence for 
the elusive T=0 pairing phase has yet been found. There 
are conflicting messages coming from calculations based 
on the quasiparticle approach. In some models, the T=0 
and T—1 pairing modes are mutually exclusive, while in 
others they are not. What is clear, however, that pre- 
dictions of calculations that impose some symmetry con- 
straints (which can rule out the presence of some pairing 
fields), should be taken with the grain of salt. 

In this work, we propose the most general nuclear en- 
ergy density functional which is quadratic in isoscalar and 



isovector densities. To this end, we discuss the isospin 
structure of the density matrices and self-consistent mean 
fields that appear in the coordinate-space HFB theory 
allowing for a microscopic description of pairing correla- 
tions in all isospin channels. The resulting expressions in- 
corporate an arbitrary mixing between protons and neu- 
trons. No particular self-consistent symmetries of the 
energy density functional have been imposed, however, 
the consequences of the time reversal and proton- neutron 
symmetry are discussed. The obtained nuclear energy 
density functional (39-41) does not have to be related 
to any given local potential. However, if the underlying 
potential is local and velocity- independent, the potential 
energy density is invariant with respect to a local gauge 
transformation. The resulting densities appear in cer- 
tain gauge-invariant combinations (47,48) which lead to 
a significant simplification of the functional. 

The self-consistent wave functions obtained by solv- 
ing the generalized HFB equations are not eigenstates of 
isospin. This is a serious drawback of the quasiparticle 
approach. To cure this problem, isospin should be re- 
stored by means of, e.g., projection techniques. While 
this can be carried out in a straightforward manner for 
energy functional that are related to a two-body po- 
tential, the restoration of spontaneously broken symme- 
tries of a general density functional poses a conceptional 
dilemma [185-188] and a serious challenge that is left for 
the future work. 
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